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1  Introduction 


Coastal  regions  represent  the  boundary  where  the  effects  of  weather  systems  off-shore  impact  activ¬ 
ities  on  land.  One  of  the  main  features  of  this  boundary  is  the  shoaling  of  waves  over  the  continental 
shelf.  As  waves  shoal,  their  amplitude  increases,  resulting  in  large  loads  on  near-shore  structures, 
etc.  This  increase  in  amplitude  can  also  lead  to  wave  breaking,  which  can  drive  near-shore  circula¬ 
tion,  with  the  resulting  currents  playing  a  significant  role  in  sediment  transport  and  beach  erosion. 
In  addition,  for  amphibious  military  operations,  the  near-shore  region  with  its  steep  and  breaking 
waves,  and  attendant  currents,  represents  a  natural  barrier  of  sorts  which  needs  to  be  traversed 
in  an  efficient  fashion.  For  all  of  these  reasons,  there  is  an  abiding  interest  in  the  monitoring  and 
prediction  of  near-shore  ocean  waves. 

This  report  describes  a  approach  which  can  be  used  to  take  time-histories  of  velocity  and  sea- 
surface  elevation  from  a  number  of  discrete  spatial  locations  and  use  them  to  produce  a  synoptic 
estimate  of  the  wave  field  for  the  region  of  interest.  The  approach  relies  on  a  variational  data 
assimilation  algorithm  which  makes  use  of  the  extended  Boussinesq  model  of  Wei  et  al  (1996).  In 
this  document,  the  mathematical  framework  for  assimilation  is  first  presented.  It  assumes  that  the 
Boussinesq  model  is  an  exact  representation  of  the  physics  of  near-shore  waves,  and  seeks  a  solution 
to  the  model  which  best  fits  the  data.  The  key  input  to  the  Boussinesq  model  which  produces  the 
desired  solution  is  the  wavemaker  source  function.  To  determine  the  source  function  which  yields 
the  best  fit  to  the  data,  a  cost  function,  a  measure  of  error  in  the  predictions,  is  first  defined.  A  set 
of  adjoint  equations  is  then  developed.  The  solution  of  the  adjoint  equations  is  used  to  calculate  the 
gradient  of  the  cost  function  with  respect  to  the  wave  maker  source  function.  This  gradient  is  used 
in  an  iterative  minimization  scheme  to  determine  the  source  function  which  yields  the  minimum 
error. 

In  this  program,  a  parallel  FORTRAN  code  was  developed  that  implements  the  assimilation 
procedure.  As  part  of  the  code  development,  an  existing  parallel  Boussinesq  code  based  on  both 
the  formulation  and  numerical  approach  of  Wei  et  al.  (1996),  was  extended  to  include  modeling  for 
run-up,  wave-breaking  and  bottom  friction.  In  addition  both  monochromatic  and  broad-spectrum 
wavemakers  were  implemented.  This  results  in  a  Boussinesq  model  code  capable  of  generating 
realistic  wave  fields  either  as  a  product  of  the  assimilation  process,  or  as  a  direct  forward  prediction 
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based  on  a  known  incident  wave  spectrum. 

The  assimilation  procedure  was  applied  to  two  cases,  observations  of  monochromatic  waves 
and  irregular  waves.  Application  of  the  procedure  for  monochromatic  waves  converged  quickly 
and  comparison  of  the  estimated  time  series  at  the  observation  locations  to  the  observations  shows 
good  agreement  between  the  predictions  and  observations.  Application  for  irregular  waves  produces 
encouraging,  but  less-accurate  results. 

In  the  following  sections,  the  assimilation  approach  is  developed,  followed  by  a  description  of 
the  Boussinesq  model  extensions  for  wave  run-up,  breaking  and  bottom  friction.  Then,  simulated 
data  to  be  used  in  the  assimilation  procedure  is  presented.  Application  of  the  procedure  to  two 
different  situations,  monochromatic  waves  and  irregular  waves,  is  then  discussed.  The  conclusions 
of  the  study  are  then  presented,  along  with  recommendations  for  future  work. 

2  The  extended  Boussinesq  model 

The  extended  Boussinesq  model  of  Wei  et  al.  (1996)  forms  the  basis  for  this  work.  This  model 
consists  of  approximate  forms  of  the  continuity  and  momentum  equations  for  water  of  finite  depth. 
The  solution  yields  the  surface  elevation  77(x,  t),  where  x  =  (x,  y)  is  position  in  the  horizontal  plane, 
and  the  horizontal  velocity  u(x, t)  =  (u,v)  at  a  reference  depth  z  =  2Q.  The  reference  depth  is 
defined  as  zQ  =  0oh  where  0$  =  -  0.531  and  h  =  /i(x)  is  the  water  depth. 

2.1  Governing  equations 

For  the  work  in  this  effort,  the  weakly  nonlinear  form  of  the  model  is  used  (after  Wei  et  al.,  1996). 
It  consists  of  a  continuity  equation 


Vt 

=  E(r},u,v)  +  S(x,f)  , 

(1) 

and  two  momentum  equations 

=  E(t),  u,  v)  +  [Fi  (t/)]t  +  F3  , 

(2) 

[V(v)]t 

=  G(r],u,v)  +  [Gi(u)]j  +  4-G5  , 

(3) 

2 


where 


u  = 

u  +  \bih2uxx  +  b2 h, (hu)xx^  , 

(4) 

V  = 

v  +  [bih2vyy  +  b2h  (hv)yyj  , 

(5) 

E  = 

-[(*  +  f?)*4  "  [(h  +  v)v]y 

-  j/i  jai h2  (uxx  +  vxy)  +  a2h  ( hu)xx  +  a2h  j 

-{h  [ai/12  ( uxy  +  Vyy)  +  a2h  (/m)xy  +  a2h  (/™)yy]  }y  , 

(6) 

F  = 

~9Vx  ~  {uux  +  vuy)  , 

(7) 

G  = 

~9Vy  ~  (uvx  +  vvy)  , 

(8) 

Fi  = 

-h  [bihvxv  +  b2  (/w)xy]  , 

(9) 

Gi  = 

h  hiiixy  b2 

(10) 

In  addition,  5(x,£)  in  (1)  is  an  internal  ‘wavemaker’  source  function,  and  Fs  and  Gs  in  (2)  and  (3), 
respectively,  are  source  terms  representing  wave- breaking,  bottom  friction,  etc.;  these  are  discussed 
below.  In  the  above  expressions,  a\  =  /?q/ 2  —  1/6,  0,2  =  Po  —  1/2,  b\  =  0\j 2  and  62  =  Po\  where 
Po  =  za/h=  -0.531. 

The  boundary  conditions  for  these  equations  vary  depending  upon  the  situation.  For  a  domain 
surrounded  by  an  impermeable,  perfectly  reflecting  boundary,  the  velocity  and  the  surface-elevation 
gradient  in  the  boundary-normal  direction  vanish  at  the  boundary.  Waves  can  be  generated  using 
a  prescribed  mass  ‘source’  in  a  region  which  acts  as  a  wavemaker.  Alternatively,  prescribed,  time- 
varying  values  of  velocity  and  surface  elevation  on  the  boundaries  can  account  waves  entering  the 
domain  from  outside. 

2.2  Additional  physical  models 

The  preceding  equations  apply  to  propagation  of  waves  in  a  region  of  finite  water  depth,  but  do  not 
adequately  treat  beaches,  where  the  water  depth  goes  to  zero  and  the  waves  run  up  on  the  beach, 
or  breaking  waves  where  steep  waves  lose  energy  through  dissipation  while  momentum  is  conserved. 
Additional  models  are  necessary  to  capture  these  physical  effects.  In  addition,  to  simulate  real- 
world  situations,  a  means  for  generating  waves  at  the  domain  boundary  is  necessary.  In  this  section, 
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wave  run-up  and  wave-breaking  models  used  in  the  Boussinesq  model  are  described,  along  with  the 
wave- generation  approach  used. 

2.2.1  Wave  run-up  model 

To  accommodate  wave  run-up  on  a  beach,  a  ‘slot7  run-up  model  is  used.  This  model  effectively 
allows  the  water  surface  to  exist  for  the  entire  computational  domain  (including  regions  of  dry 
beach).  This  is  accomplished  by  relegating  the  water  to  narrow  slots  (or  equivalently,  porous 
regions)  in  the  ‘dry’  portions  of  the  domain.  In  the  dry  regions,  the  amount  of  water  is  kept 
negligibly  small,  but  non-zero.  As  a  wave  crest  approaches  the  beach,  the  slot  rapidly  fills  and  the 
water  level  rises  above  the  bottom  level;  as  the  wave  recedes,  the  water  runs  off  until  only  the  small 
amount  that  resides  in  the  slot  remains.  The  basic  approach  is  described  by  Kennedy  et  al.  (2000) 
and  Chen  et  al.  (2000);  a  slightly  modified  approach  is  used  here. 

The  slot  run-up  method  is  implemented  using  a  modified  continuity  equation,  based  on  the 
exact  equation  (1),  above.  The  modified  equation  is 

Vt  =  -  {Mv)u  +  h  [ai h2  (uxx  +  vxy)  +  a,2h  ( hu)xx  +  a2h  (ftv)^]  j 

-  |.4( T])v  +  h  [ai/i2  (Uxy  +  Vyy)  +  a2h  {hu)^  +  a2h  {hv)yyj  |  +  S(x,  t )  ,  (11) 

where  the  run-up  model  is  represented  by  A(rj ).  For  no  run-up  modeling,  A(tj)  —  h  +  77;  the  model 
is  invoked  by  setting 


A(rj)  =  hc  +  V  (/i  +  T|)  <  0 

=  (ft  T  77)  4*  {hc  —  h)  8  (h  +  7j)>0,  (12) 

where  6  is  the  relative  cross-sectional  area  of  slots  ((5  is  typically  set  to  0.02),  and  hc  is  set  to  the 
maximum  water  depth  in  the  domain. 

This  implementation  of  the  slot  rumup  model  undergoes  an  abrupt  transition  as  the  waves  run 
up  and  ‘fill’  the  slot,  compared  to  the  form  of  Kennedy  et  al.  (2000),  which  has  a  gradual  transition; 
however,  in  practice  a  gradual  transition  does  not  ensure  smoothness  in  the  horizontal  derivatives 
of  A (77),  and  lack  of  smoothness  can  lead  to  instability.  For  this  reason,  a  high- limit  filter  which 
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cuts  off  at  twice  the  grid  spacing  is  applied  to  A{rj),  once  it  is  calculated.  This  ensures  well  behaved 
derivatives  and  results  in  a  robust  model. 

2.2.2  Wave  breaking  and  bottom  friction 

Wave  breaking  is  modeled  in  a  manner  similar  to  that  described  in  Kennedy  et  al.  (2000),  i.e.  by 
adding  an  additional  diffusive  term  to  the  momentum  equations.  This  is  given  by 

Fs,b  =  +  +  2(/t  I|_  q)  fr  +  +  Ub  ^  +  (13) 

G°’b  =  p6  [(ft  +  n)v\y}y  +  2(fe1+  fr  [(h  +  v)v)x  +  vb  [(ft  +  r1)u\y)x  (14) 

where  Vf,  is  the  eddy  viscosity  associated  with  breaking.  This  eddy  viscosity  is  defined  after  Zelt 
(1991)  as 

Vb  =  B5l{h  +  r))r)t  (15) 


where  St,  =  1.2  and 


B  =  min 

max  (o,  %  -  l)  ,  1 

V  Vt  )  \ 

(16) 


where  ?7£*  =  0.3 \/g(h  +  rf).  The  form  of  (16)  causes  B  to  be  bounded  by  zero  and  one  and  vary 
linearly  between  those  limits  as  rjt  increases  from  rf  to  2rf. 

The  effects  of  bottom  friction  are  modeled  via  source  terms  in  the  momentum  equations  which 
oppose  the  fluid  motion.  The  model  form  is  given  by 


II 

r,  luF 
h  +  ij 

(17) 

GsJ  — 

~  lu|v 

h  +  r) 

(18) 

where  Q,  =  0.005. 
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2.2.3  Wavemaker  source  function 


Monochromatic  waves  are  generated  using  a  source  term  in  the  continuity  equation  (1),  as  described 
in  Wei  et  al.  (1999).  The  source  function  amplitude  is  set  to 


UJ2  -  (oo  +  %)gk(kho)z 
uak  [1  —  ao(A;/io)2] 


2?7o  cos($) 


(19) 


where  770  is  the  wave  amplitude,  u)  is  the  wave  circular  frequency,  9  is  the  wave  direction  measured 
counter-clockwise  from  the  x  axis,  k  =  2ir/\  is  the  wavenumber,  ho  is  the  fixed  depth  under  the 
wavemaker,  and  a  =  b\  +  &2-  Here,  a  is  given  by 


a  = 


(20) 


where  kx  A;cos(0)  and  j3\  —  5(A/2)  2.  The  source  function  used  in  (1)  for  a  wavemaker  located 
at  x  —  0  on  the  left-hand  boundary  of  the  computational  domain  is  then  given  by 


S(xyt)  =  w(x)S{y1t) 

_  e~P*x2  Sq  sin (kyy  —  ut  T  <j>) 


(21) 


where  tu(a;)  =  e~®l'x2 ,  ky  =  ksin(0)  and  0  <  4>  <  2n  is  an  assigned  phase. 

For  irregular  waves,  the  spectral  densities  of  the  Fourier  modes  comprising  the  wave  spectrum 
are  converted  to  equivalent  amplitudes  and  treated  as  monochromatic  waves.  The  wavemaker 
source  function  is  then  given  by 

/vm 

S{x,t)  =  Y^e~0lX  ■sosin(A:’y  -  u/f  +  4>l)  ,  (22) 

i 

where  Nm  is  the  number  of  modes  and 


(w*)2  -  (qq  +  Dfffc^fc^o)3 
w%alkl  [1  —  a(,{k'Lho)2\ 


2?7o  cos(9l) 


(23) 
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and 


-Hf  /m  '  (24) 

Here,  the  value  of  (3\  is  calculated  from  the  dominant  wave  (spectral  peak),  resulting  in  a  fixed 
wavemaker  width  for  all  spectral  components. 

The  present  implementation  allows  the  specification  of  the  wave  amplitude,  period  and  direction 
from  the  command  line  for  monochromatic  waves.  For  irregular  waves,  a  frequency-direction 
spectrum  file  in  the  FRF  8m  Array  format  is  read. 

2.3  Numerical  implementation 

The  Boussinesq  equations  incorporating  the  physical  models  described  above  are  solved  using  a 
parallel  FORTRAN  code.  The  numerical  implementation  is  identical  to  that  described  in  Wei 
et  al.  (1996):  explicit  time  integration  using  a  third  order  predictor  and  an  iterated  fourth  order 
corrector.  The  spatial  grid  is  Cartesian  with  uniform  spacing,  and  the  grid  is  generated  internally 
by  the  code.  The  solution  is  filtered  every  twenty  time  steps  to  eliminate  spurious  short-wavelength 
components  (which  result  from  numerical  errors  and  nonlinearities  and  can  accumulate  due  to  the 
absence  of  natural  diffusion  or  dissipation  in  the  equations).  The  code  is  parallelized  using  domain 
decomposition  and  the  message-passing  interface  (MPI).  It  will  run  for  an  arbitrary  number  of 
processors  (the  number  of  processors  is  specified  on  the  command  line)  on  both  Linux  (Portland 
Group  compilers  and  LAM  MPI)  and  SUN  Microsystems  (Forte  compiler  suite  and  MPICH  MPI) 
platforms.  Bathymetry  information  is  read  from  a  binary  file  on  the  same  grid  as  the  computational 
grid,  while  the  tide  is  set  via  a  command  line  argument. 

3  The  assimilation  methodology 

It  is  desired  to  determine  the  wavemaker  source  function  which  will  yield  predictions  from  the 
extended  Boussinesq  model  which  match  most  closely  a  set  of  time-series  observations  of  the  near¬ 
shore  surface  elevation  and  velocity.  The  approach  taken  will  be  based  on  the  variational  method 
described  by  Le  Dimet  &  Talagrand  (1986),  corresponding  to  the  ‘strong-constraint’  formalism 
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described  by  (Bennett,  1992,  p.  148).  The  objective  is  to  determine  the  source  function  which 
minimizes  a  ‘cost  function’  (a  positive-definite  measure  of  the  error  in  the  predictions)  subject 
to  the  constraint  the  the  dynamical  equations  governing  the  wave  field  (the  extended  Boussinesq 
equations)  are  exactly  satisfied. 

The  above  constrained  minimization  problem  is  equivalent  to  the  unconstrained  minimization 
of  the  cost  function  ‘augmented’  with  the  product  of  the  constraint  equations  and  a  set  of  ‘adjoint’ 
variables  (i.e.  Lagrange  multipliers).  The  cost  function  will  be  minimized  for  the  set  of  model 
inputs  which  correspond  to  a  stationary  point  of  this  augmented  cost  function;  i.e.  where  its  first 
variation  vanishes.  The  requirement  that  the  first  variation  vanish  yields  the  governing  equations 
for  the  adjoint  variables,  as  well  as  the  relationship  between  the  adjoint  variables  and  the  gradient 
of  the  cost  function  with  respect  to  the  control  variables,  i.e.  the  wavemaker  source  function. 

In  the  following  sections,  the  cost  function  and  the  augmented  cost  function  are  described. 
The  first  variation  of  the  cost  function  is  taken  and  the  conditions  under  which  the  first  variation 
vanishes,  the  Euler-Lagrange  equations,  are  derived.  The  Euler-Lagrange  equations  are  then  used 
to  develop  both  the  adjoint  equations  and  the  relationship  between  the  adjoint  solution  and  the 
gradient  of  the  cost  function  with  respect  to  the  wavemaker  source  function. 

3.1  The  cost  function 

The  cost  function  is  comprised  of  two  parts.  The  first  part  is  a  measure  of  the  error  in  the  predicted 
wave  field.  The  second  part  is  based  on  the  model  input  being  estimated  (i.e.  the  wavemaker  source 
function),  and  is  included  to  ensure  a  unique  solution. 

The  first  part  of  the  cost  function  is  a  positive-definite  measure  of  the  error  in  the  predicted 
wave  field.  Observations  are  taken  to  consist  of  time  records  of  length  T  of  the  surface  elevation 
and  vector  velocity  at  N  discrete  locations  in  the  domain.  In  equation  form,  this  is  given  by 

J\  =  E/  /  kl  {v  ~  Vi)2  +  V2—  (u  -  2»)2  +</?3~  (v  -Vif 
i  Jo  JnY  9  9 

where  the  circumflex  indicates  an  observed  quantity,  the  delta  functions  specify  the  locations  of  the 
observations,  and  the  (pk  are  constants  which  allow  variable  weighting  of  the  different  observation 
quantities.  The  integration  here  is  over  the  entire  space-time  region  of  interest,  0  <  t  <  T  and  for 


<>(X  -  Xj) 

NT 


dxdt  ,  (25) 
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all  spatial  locations  x  €  71. 

The  second  part  of  the  cost  function  is  based  on  the  estimated  model  inputs,  in  this  case  the 
wavemaker  source  function.  In  order  to  assure  a  unique  solution,  the  input  variable(s)  must  also 
be  penalized  via  inclusion  in  the  cost  function  (Bennett  &  Miller,  1991).  This  can  be  accomplished 
by  including  one  or  more  positive-definite  measures  of  the  wavemaker  source  function,  which  leads 
to 


h 


-  rj 

Jo  Jb 


h  o 


b  gLsT  _ 


d£dt 


(26) 


where  £  indicates  position  along  the  centerline  of  the  wavemaker  B  (which  could  correspond  to  a 
domain  boundary),  Ls  is  the  length  of  the  wavemaker,  and  Sc  is  the  slope  of  the  source  function 
along  the  wavemaker.  The  first  term  constrains  the  source  function  to  have  minimum  possible 
energy,  while  the  second  constrains  the  directional  spread  of  the  generated  wave  spectrum. 

Combining  J\  and  J2  yields  the  complete  cost  function: 


J  = 


J0  fn\v  Vi)2  +  ¥>2^7  («  -  2i)2  +  <p  2^  ( v  -  Vi)2 


+ 


in 

rT  r  ho 


5(x  -  Xi) 
NT 


dxdt 
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Jo  Jb 


b  gLsT 


<P3 S  {t,t)+<p4hQSfa,t) 


d£dt 


(27) 


The  methods  of  variational  calculus  (see  e.g.  Courant  k,  Hilbert,  1953)  show  that  minimizing  the 
cost  function  J  subject  to  the  constraint  that  the  solution  satisfy  the  extended  Boussinesq  equations 
is  equivalent  to  unconstrained  minimization  of  an  ‘augmented’  cost  function.  The  physical  solution 
variables  are  (77,  u ,  v)  and  the  corresponding  governing  equations  axe  given  by 


({77,  -  E  -  wS}  ,  {Ut  -  F  -  Fit]  ,  {Vt  -  G  -  Gu})  =  0.  (28) 

(Here,  the  source  terms  in  the  momentum  equations  are  ignored  in  the  development  of  the  adjoint 
equations.)  A  set  of  adjoint  variables  (&,/?,  7)  is  introduced  which  correspond  to  (r)tu,v).  The  cost 
function  is  augmented  by  the  product  of  these  adjoint  variables  and  the  governing  equations  (the 
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extended  Boussinesq  model).  This  augmented  cost  function  £  is  given  by 


£ 


-E-wS}+P{Ut-F-  Flt}  +7  {Vt-G-  Gu}]dxdt  . 


(29) 


Here,  £  is  a  functional  of  a,  P,  7,  77,  u)  v  and  S,  i.e.  £  =  £(a,/3,7,77,  u,v,  S). 


3.2  The  Euler-Lagrange  equations 


The  Euler-Lagrange  equations  define  the  conditions  for  the  minimum  of  the  augmented  cost  func¬ 


tion.  The  cost  function  will  be  at  a  minimum  for  a  stationary  point,  i.e.  where  the  first  variation 


of  £  the  vanishes.  The  first  variation  is  given  by 


C/1  dC  J  3£  dC  J  dC  dC  J  dC  J  dC 
SC  =  —da  +  —dp  +  — d7  +  —  drj  +  —  du  4-  ^-dv  +  —dS 
oa  dp  oj  or]  on  ov  oS 


(30) 


and  SC  —  0  for  a  stationary  point.  Hence,  each  of  the  terms  in  (30)  must  vanish  independently, 


which  requires  that  the  partial  derivatives  of  £  all  be  zero:  The  resulting  set  of  requirements  that 


this  imposes  are  the  Euler-Lagrange  equations.  These  are  examined  in  the  next  sections. 


3.2.1  The  constraint  equations 

Setting  the  first  three  terms  in  (30)  to  zero  yields  the  condition  that  the  constraint  equations,  i.e. 
the  extended  Boussinesq  model,  must  be  satisfied  at  the  minimum  in  £.  These  terms  are 


—da  = 
da 

f  [  [{rjt  —  E  —  wS }  Sa\  dxdt  , 

Jo  Jn 

(31) 

>  - 

[  [  [{Ut-F-Fu}6p]dxdt, 

Jo  Jn 

(32) 

-K~dl  = 
<77 

[  [  \{Vt-G-Gu}5'y]dxdt. 

Jo  Jn 

(33) 

Since  the  regions  for  the  integrations  are  arbitrary,  the  integrands  must  be  identically  zero  for  the 
integrals  to  vanish.  In  addition,  the  values  of  (5a,  Sp  and  <57  are  arbitrary,  the  terms  in  braces 
must  vanish  independently.  This  requirement  is  the  same  as  requiring  the  extended  Boussinesq 
equations  (1-3)  to  be  satisfied  exactly.  (This  is  the  ‘strong’  constraint  on  the  wave  field.)  This 
essentially  requires  that  the  wave  field  which  gives  the  minimum  cost  function  (agrees  best  with 
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the  observations)  has  the  appropriate  dynamics. 


3.2.2  The  adjoint  Boussinesq  equations 

The  next  three  terms  of  (30),  those  involving  Sr],  Su  and  Sv ,  yield  a  set  of  equations  for  the  adjoint 
variables  a,  (3  and  7.  The  first  variation  of  £  with  respect  to  77  is  given  by 


5r) 


dxdt  , 


(34) 


Here,  the  U>  V ,  F\  and  G 1  terms  do  not  appear  since  they  have  no  rj  dependence.  After  integrating 
by  parts  to  move  the  Srfs  outside  the  derivatives  and  rearranging,  the  result  is 


at~E'  ~  NT  _  ^ 


Sr]  dxdt  , 


(35) 


where  the  final  term  in  the  brackets  comes  from  the  77-dependent  term  in  the  cost  function.  The 
first  variation  of  (30)  will  vanish  independent  of  the  value  of  Sr]  only  if  the  term  in  brackets  is  zero. 
This  yields  the  governing  equation  for  a,  the  adjoint  of  the  variable  rj: 


where 


at  =  E>  +  NT  “  Xi 


(36) 


E*  =  -uax  -  vay  -  g(/3x  +  7y)  . 


(37) 


To  obtain  the  equation  for  (3,  the  adjoint  of  tz,  we  take  the  first  variation  with  respect  to  u . 
This  gives  the  result 


8C.  dJ  J  fTf  [  f  dE\  _  a(dUt  dF  \,  /  8G  dGlt 


Su 


dx  dt  .  (38) 


After  integrating  by  parts  and  collecting  terms  this  yields 


=  -  JoJn  [U'iP) ]|  -  F'  -  [T^lt  -  XA  -  ^)<5(x  -  x*) 


-5- du 

ou  JoJn 


Su  dx  dt 


(39) 
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The  first  variation  of  (30)  will  vanish  independent  of  the  value  of  Su  only  if  the  term  in  brackets  is 
zero.  This  yields  the  governing  equation  for  /?,  the  adjoint  of  the  variable  u.  The  equation  is 

n  =  (40) 

i 

where 


u'  =  p+[b,{h2p) 

XX  +  b2h(hp)  XX  ]  (41) 

F  —  +  n)ax  T  CLi  0-x)xx  4~  ( h  T  ^ ®x)xx  T  {h 

+uxf3  -  ( Up)x  -  ( V0)y  4-  Vxi  (42) 

F[  =  -  [bi(h2^f)xy b2h(hj)xy^  .  (43) 

The  equation  governing  7,  the  adjoint  of  t>,  can  be  derived  by  inspection  from  the  (3  equation. 
Observing  that  the  u  and  v  equations  are  symmetric  with  respect  to  (tqu)  and  (x,y),  i.e.  inter¬ 
changing  x  and  y,  as  well  as  u  and  v}  in  the  u  equation  (2)  yields  the  v  equation  (3),  and  vice  versa. 
A  similar  interchange  in  the  rj  equation  (1),  yields  (1).  Hence,  the  governing  equations  for  7  is 

\y'\t  =  G'  +  [£i](  +  T(v  ~  ?i)<5(x  -  xi)  (44) 

i 

where 

V  =  7  +  [M/i2  l)yy  +  b2h{h-y)yy\ 

Cr  —  j^(/l  4“  TpCty  “f"  fll  |(/l  ^y)yy  4”  (^  ^x)xj/)  4"  Oj2h  ^(/i  Oiy)yy  ~\~  (/l  ^x)jylj 

4-^7  -  {^l)v  -  («7)x  4-  UyP 
G\  =  -  [bi(h2P)xy  +  b2h(hP)Xy]  . 

3.2.3  The  cost  function  gradient 

This  final  term  in  (30)  involves  the  wavemaker  source  function  S(£,t),  where  £  is  the  coordinate 
along  the  length  of  the  wavemaker.  (Here,  x  =  (£,  £),  where  £  is  the  coordinate  along  the  wavemaker 
and  £  is  the  coordinate  normal  to  the  wavemaker;  for  the  work  presented  here,  (C,£)  =  (rr,  2/) T  but 


(45) 

(46) 

(47) 
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this  is  not  general  requirement.)  If  the  constraint  (Boussinesq)  equations  and  the  adjoint  equations 
are  satisfied,  then  last  term  in  (30)  is  the  only  one  that  remains,  and  therefore 


■ 


(48) 


This  quantity  2/§dS  contains  contributions  from  the  cost  function  and  the  continuity  equation, 
where  the  wavemaker  source  S  appears: 


5C 


dJ 

dS 

rT 


dS  + 


wSS 


dx  dt 


/  7 

Vo  Jb 


B  9L3T 


ai-° 

S  SS  ^4^0^^  SS^ 


2ho 


a  w(C)  SS 


dx  dt 


(49) 


Integrating  by  parts  and  collecting  terms  yields 


SC  = 


=  JQ  jB{^  -^a(C,0  rn(C)dc}  6S  d$dt  , 


(50) 


where  the  integral  inside  the  braces  is  over  lines  normal  to  the  wavemaker  centerline  (indicated  by 
AT),  i.e.  in  the  £-direction. 

Equation  (50)  indicates  that  a  small  decrease  in  C  (the  thing  we  are  trying  to  minimize)  will 
result  from  a  change  in  the  source  function  SS  given  by 


SS  oc  — 


f  2h0 
\gLsT 


ipzS  -  (p4hl S# 


jf_a(C,0  tu(C)dc)  • 


(51) 


Hence,  the  gradient  of  the  augmented  cost  function  C  with  respect  to  the  wavemaker  source  function 
S  is  effectively  given  by 


dC 

dS 


2  h0 
gLsT 


S  —  (p4hoStf 


[  <*( CiO  MCMC  • 

JAf 


(52) 


This  gradient  can  be  used  to  determine  the  value  of  the  wavemaker  source  function  which  will  yield 
the  best  agreement  between  the  estimated  wave  field  and  the  observation  data. 

Examination  of  (52)  shows  that  the  gradient  has  three  components.  The  first  two  depend  on  the 
source  function  S  and  its  second  derivative  %.  and  are  included  to  produce  a  unique  solution.  The 
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third  depends  on  the  adjoint  solution  integrated  over  the  region  of  the  wavemaker.  This  portion  of 
the  gradient  relates  directly  to  the  error  in  the  predicted  wave  field  . 

3.3  Solution  of  the  adjoint  equations 

The  adjoint  equations  are  similar  in  mathematical  character  to  the  extended  Boussinesq  equations 
(they  also  have  identical  linear  dispersion  properties).  As  a  result,  the  same  solution  techniques 
that  are  used  to  solve  the  extended  Boussinesq  model  equations  can  be  applied  to  solve  the  adjoint 
equations.  While  the  boundary  and  initial  conditions  for  the  adjoint  equations  are  somewhat 
problem-dependent,  a  few  general  statements  can  be  made. 

The  adjoint  equations  are  solved  for  the  time  period  for  which  the  observations  are  available. 
They  are  integrated  backward  in  time  from  the  final  time,  t  =  T.  The  ‘initial’  conditions  for  the 
adjoint  are  ot(x,T)  =  /?(x,T)  =  7 (x,  T)  =  0. 

The  spatial  boundary  conditions  for  the  adjoint  equations  are  taken  to  be  the  same  as  for 
forward  equations:  a  has  zero-gradient  conditions  at  all  boundaries;  (3  =  0  and  7  =  0  a  the  x  and 
y  boundaries,  respectively,  and  zero-gradient  conditions  elsewhere. 

Each  adjoint  equation  contains  an  ‘observation’  term  for  input  of  observations  of  the  correspond¬ 
ing  physical  variable  (77  for  the  a  equation,  u  for  the  (3  equations,  and  so  on).  Since  the  adjoint 
equations,  and  their  boundary  and  initial  conditions  are  otherwise  homogenous,  the  observation 
terms  are  what  make  the  adjoint  solution  non-zero.  When  the  predictions  for  77,  etc.  agree  exactly 
with  the  observations,  the  observation  terms  become  zero,  and  the  resulting  solution  to  the  adjoint 
equations  is  identically  zero.  The  gradient  of  the  cost  function  with  respect  to  the  inputs  for  the 
Boussinesq  model  will  then  be  zero,  and  the  assimilation  procedure  will  have  converged. 

3.4  Assimilation  algorithm 

The  cost  function  gradient  calculated  from  the  adjoint  solution  is  used  in  a  conjugate-gradient 
minimization  scheme  where  cost  function  is  minimized  by  adjusting  the  wavemaker  source  function. 
The  procedure  goes  as  follows: 

1.  The  adjoint  equations  are  solved  backward  in  time  with  the  observations  as  input  and  the 
gradient  of  the  cost  function  with  respect  to  the  wavemaker  source  function  is  calculated  from 
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the  adjoint  solution  in  the  region  of  the  wavemaker  using  (52). 

2.  The  direction  of  steepest  descent  (opposite  the  gradient)  is  defined  as  the  search  direction. 

3.  The  location  of  the  minimum  for  the  cost  function  along  the  search  direction  is  determined  us¬ 
ing  a  secant  method  to  find  the  location  where  the  directional  derivative  (the  inner  product  of 
the  local  gradient  with  the  search  direction)  vanishes.  Each  iteration  requires  a  forward  model 
solution  with  the  updated  source  function,  and  an  adjoint  model  solution  with  the  forward 
solution  as  input,  followed  by  the  calculation  of  the  gradient  from  the  adjoint  solution. 

4.  The  most  recent  gradient  and  the  original  search  direction  are  used  in  the  determination  of 
the  new  search  direction  using  the  conjugate  gradient  scheme. 

The  last  two  steps  represent  a  conjugate  gradient  cycle,  and  are  repeated  until  the  forward  predic¬ 
tion  converges.  Once  the  procedure  has  converged,  the  most  recent  forward  prediction  represents 
the  best  fit  to  the  observation  data.  The  final  wavemaker  source  function  is  saved  to  a  file  and 
can  be  used  as  input  to  the  forward  model  to  reproduce  the  best-fit  wave  prediction  outside  of  the 
assimilation  procedure. 

This  procedure  is  implemented  in  a  parallel  FORTRAN  code  which  executes  both  the  extended 
Boussinesq  model  and  its  adjoint  and  controls  the  conjugate  gradient  minimization  algorithm.  The 
code  is  parallelized  using  domain  decomposition  and  the  message- passing  interface  (MPI).  As  in  the 
description  of  the  Boussinesq  model  above,  it  will  run  for  an  arbitrary  number  of  processors  (the 
number  of  processors  is  specified  on  the  command  line)  on  both  Linux  (Portland  Group  compilers 
and  LAM  MPI)  and  SUN  Microsystems  (Forte  compiler  suite  and  MPICH  MPI)  platforms.  Inputs 
are  the  same  as  for  the  Boussinesq  model  are  described  above  with  the  addition  of  a  data  file  for 
each  of  the  data  records  used  in  the  assimilation.  Execution  of  the  forward  model  or  the  assimilation 
procedure  is  controlled  via  command-line  arguments,  as  are  the  parameters  of  the  iteration  scheme. 

4  Forward  Boussinesq  model  calculations 

The  extended  Boussinesq  model  code  including  the  wave-breaking  and  run-up  models  and  the 
wavemaker  source  function  was  tested  by  calculating  three  different  situations.  The  first  was 
monochromatic  waves  propagating  across  a  realistic  bathymetry,  derived  from  a  survey  for  a  1  km 
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square  area  around  the  US  Army  Corps  of  Engineers  Field  Research  Facility  (FRF)  in  Duck,  NC. 
This  tested  the  basic  wave  propagation  capabilities  of  the  code  and  served  as  validation  of  the  run¬ 
up  model  in  particular  (a  run-up  implementation  that  works  on  a  planar  beach  is  not  necessarily 
robust  enough  for  a  realistic  beach  topography).  The  second  case  was  the  same  bathymetry,  but 
for  irregular  waves  generated  from  a  measured  FRF  wave  spectrum.  This  tested  the  model  for 
propagation  and  run-up  of  multiple  wave  modes,  and  also  verified  the  broad  spectrum  directional 
wavemaker  implementation. 

4.1  Monochromatic  waves 

To  test  the  basic  wave  generation,  propagation  and  run-up  in  the  Boussinesq  model  code,  it  was 
applied  to  a  1  km  square  region  around  the  USACE  FRF  in  Duck,  NC.  The  model  was  run  on  a 
2  m  x  10  m  (cross-shore  by  along-shore)  grid.  The  bathymetry  is  shown  in  figure  1,  and  is  seen  to 
vary  from  zero  to  8  m,  with  a  bar  located  roughly  200  m  offshore.  There  is  a  break  in  the  bar  at  the 
location  the  FRF  pier.  To  test  the  assimilation  algorithm,  synthetic  measurements  were  extracted 
at  the  locations  indicated  in  the  figure,  roughly  at  the  crest  of  the  bar. 

The  Boussinesq  model  was  run  for  this  bathymetry  using  a  time  step  of  0.1  s.  Waves  were 
generated  by  a  mass  source  at  the  left-hand  boundary.  The  wavemaker  source  region  centered  at 
the  boundary  ( x  =  0),  had  a  Gaussian  profile  and  extended  about  50  m  into  the  computational 
domain.  Waves  were  generated  with  a  significant  wave  height  Hmo  =  1  m,  a  wave  period  of  Tp  —  6  s 
and  by  using  a  traveling  wave  (propagating  in  the  +y-direction)  in  the  source  function  were  set  to 
propagate  at  10  degrees  to  the  cross-shore  direction.  The  resulting  wave  field  is  shown  in  figure  2. 
Here  we  can  see  that  the  waves  run  up  on  the  beach  and  the  water  surface  does  not  extend  to 
the  right-hand  side  of  the  domain.  The  wave  field  shown  in  the  figure  is  a  snapshot  of  the  surface 
elevation  after  the  waves  have  been  running  up  on  the  beach  several  minutes  of  real  time.  For 
these  calculations,  wave  absorbing  layers  are  located  within  50  m  of  the  y  boundaries  to  suppress 
reflections;  this  is  responsible  for  the  change  in  wave  direction  very  near  the  boundaries. 

In  figure  2,  seven  locations  spaced  100  m  apart  are  indicated  for  which  records  of  Tj(t ),  ti(£)  and 
v(t)  are  saved.  The  velocity  saved  is  for  the  reference  depth  2a,  about  midway  in  the  water  column, 
although  in  shallow  water  the  velocity  is  relatively  uniform  over  depth.  The  temporal  length  of 
the  records  was  several  times  the  time  required  for  a  wave  to  traverse  the  computational  domain, 
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FIGURE  1.  Bathymetry  used  in  simulations  for  a  1  km  square  area  around  the  USACE  Field 
Research  Facility  at  Duck  NC. 
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FIGURE  2.  Monochromatic  waves  with  Hmo  ~  1  m,  Tp  =  6  s  propagating  at  10  degrees  to 
the  cross-shore  direction  over  the  bathymetry  of  figure  1,  above.  The  waves  are 
generated  at  the  left,  offshore  boundary  and  run  up  on  the  beach  at  the  right. 
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FIGURE  3.  Sample  wave  records  showing  rj(t),  u(t)  and  v(t)  for  locations  3,  5  and  7  in 
figure  1,  above. 

and  the  sampling  rate  was  0.1  sec,  the  computational  time  step.  A  portion  of  three  of  those  time 
records  is  shown  in  figure  3. 

4.2  Irregular  waves 

To  calculate  realistic  wave  fields  using  the  extended  Boussinesq  model  code,  broad-spectrum  waves 
must  be  generated.  A  broad-spectrum  wavemaker,  described  above  in  section  2.2.3  was  imple¬ 
mented.  This  wavemaker  takes  as  input  a  spectrum  in  the  format  of  the  FRF  8m  array,  ASCII 
output  file.  Figure  4  shows  a  snapshot  of  the  wave  field  generated  using  a  measured  FRF  spectrum 
as  input.  The  wave  field  has  a  significant  wave  height  of  1.3  m,  a  peak  wave  period  of  5.5  s  and 
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Figure  4.  Irregular  waves  with  Hmo  =  1.3  m,  Tp  =  5.5  s  with  a  dominant  wave  direction 
roughly  cross-shore  propagating  over  the  bathymetry  of  figure  1,  above.  The 
waves  are  generated  at  the  left,  offshore  boundary  and  run  up  on  the  beach  at 
the  right. 

the  dominant  wave  direction  is  in  the  cross-shore  direction.  Figure  5  shows  a  30-minute  average  of 
the  surface  elevation  and  velocity  for  the  wave  field  shown  in  figure  4.  Time  histories  of  rj,  u  and  v 
from  locations  2,  4  and  6  in  figure  1  are  shown  in  figure  6. 

Broad-spectrum  waves  represent  a  challenge  for  nonlinear  Boussinesq  models.  The  long  waves 
with  higher  propagation  speeds  set  the  requirements  for  the  time  step  (based  on  the  local  Courant 
number).  The  short  waves  set  the  spatial  grid  requirements.  In  addition,  the  presence  of  multiple 
Fourier  components  and  the  nonlinearity  of  the  model  leads  to  the  generation  of  sum  and  difference 
frequencies.  These  can  lead  to  accumulation  of  spurious  wave  energy  in  wavelengths  near  the  grid 
spacing  (sawtooth  waves).  The  results  shown  in  figure  4,  which  represents  the  wave  field  at  the  end 
of  30  minutes  of  real-time  evolution,  shows  no  evidence  of  any  short-wave  noise. 
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FIGURE  5.  Results  from  Boussinesq  calculations  of  irregular  waves  for  a  30  min  period: 

mean  wave  height  overlaid  with  mean  velocity  vectors  (the  maximum  mean 
velocity  is  0.35  m/s)  . 
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FIGURE  6.  Time  records  of  r](t),  u(t)  and  v(t)  for  locations  2,  4  and  6  in  figure  1. 
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FIGURE  7.  The  adjoint  solution  a(x,  £)  for  the  seven  in-situ  observations  of  the  monochro¬ 
matic  wave  field.  The  gradient  of  the  cost  function  with  respect  the  wave  maker 
source  function  is  calculated  from  the  adjoint  solution  near  the  left  boundary. 

5  Assimilation  results 

In  this  section,  the  assimilation  algorithm  is  applied  using  observations  generated  using  both  the 
monochromatic  and  irregular  wave  simulations  described  above. 

5.1  Results  for  monochromatic  waves 

Time  records  of  surface  elevation  and  velocity  for  monochromatic  waves  were  assimilated  into  the 
Boussinesq  model  using  the  procedure  outlined  above.  This  was  done  for  the  wave  field  shown 
above  in  section  4.1.  The  first  step  in  the  procedure  is  to  solve  the  adjoint  Boussinesq  equations 
using  the  observations  as  input.  A  snapshot  of  the  resulting  adjoint  solution  is  shown  in  figure  7. 
This  figure  shows  the  adjoint  surface  elevation  field  a(x).  The  adjoint  solution  near  the  left-hand 
boundary  is  used  to  calculate  the  gradient  of  the  cost  function  with  respect  to  the  wavemaker  source 
function. 

Figure  8  shows  the  iteration  history  of  the  cost  function.  The  complete  cost  function  Jrotal 
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FIGURE  8.  Iteration  history  of  the  cost  function  Jrotai  (27)  for  monochromatic  waves  along 
with  the  contributions  from  the  observations  J0&s  (25)  and  the  source  function 
Ja  (26);  here  each  iteration  corresponds  to  a  conjugate-gradient  cycle.  This 
results  shows  that  the  cost  function  is  reduced  substantially  and  converges  in 
less  that  ten  iterations. 

defiend  by  (27),  above,  is  shown,  along  with  the  components  attributed  to  the  observations  J0f,s 
(25)  and  that  attributed  to  the  wavemaker  source  function  J3  (26).  Each  ‘iteration’  in  the  figure 
shows  the  multiple  values  that  are  obtained  during  a  conjugate-gradient  cycle.  As  the  iterations 
proceed,  the  portion  of  the  cost  function  attributed  to  the  observations  decreases,  while  that  from 
the  source  function  increases,  as  is  expected.  The  overall  cost  function  J  decreases  and  levels  out 
after  about  5  iterations,  and  is  clearly  converged  by  ten  iterations. 

The  estimated  sea  surface  resulting  from  the  assimilation  procedure  is  shown  in  figure  9.  With 
the  exception  of  some  small  undulations  in  the  wave  crests,  the  wave  field  is  reproduced  with 
good  accuracy.  Figure  10  shows  a  comparison  between  the  time-series  observations  used  in  the 
assimilation  process  and  those  from  the  resulting  estimate.  Again  the  agreement  is  quite  good. 
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Figure  9.  Sea  surface  elevation  estimated  using  the  assimilation  procedure;  this  result 

compares  well  to  the  actual  sea  surface  shown  in  figure  2 . 

5.2  Results  for  irregular  waves 

Time  records  of  surface  elevation  for  the  irregular  waves  of  section  4.2  were  assimilated  into  the 
Boussinesq  model  using  the  procedure  outlined  above.  The  iteration  history  of  the  cost  function 
is  shown  in  figure  11.  Here  it  can  be  seen  that  the  cost  function  reduction  is  less  than  for  the 
monochromatic  waves,  shown  above,  and  there  is  only  about  a  factor  of  two  reduction. 

A  snapshot  the  estimated  wave  field  is  shown  in  figure  12,  compared  to  the  corresponding 
snapshot  from  the  actual  wave  field.  Comparison  of  the  two  images  shows  that  while  the  long-wave 
components  of  the  sea  surface  are  captured  well,  the  short-wave  components  are  missing  from  the 
estimated  sea  surface.  This  is  particularly  true  of  short-wave  components  propagating  in  the  along¬ 
shore  direction.  Figure  13  shows  a  comparison  between  the  time-series  observations  used  in  the 
assimilation  process  and  those  from  the  resulting  estimate.  Overall  the  lower  frequency  component 
are  captured,  but  the  higher  frequency  components,  in  particular  in  the  v  velocity  are  missed.  The 
level  of  error  in  the  u  velocity  estimate  is  comparable  to  that  in  the  v  velocity  indicating  that 
oblique  waves  are  not  being  well  captured. 
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Figure  10.  Comparison  of  estimated  and  actual  time  histories  at  the  observation  locations. 


FIGURE  11.  Iteration  history  for  the  cost  function  Jrotai  (27)  along  with  the  contributions 
from  the  observations  J0^s  (25)  and  the  source  function  J3  (26);  here  each 
iteration  corresponds  to  a  conjugate-gradient  cycle. 
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Figure  12.  a)  Sea  surface  elevation  estimated  using  the  assimilation  procedure,  b)  Actual 
sea  surface  corresponding  to  the  time  for  the  estimated  sea  surface.  The  es¬ 
timated  wave  field  captures  the  longer  wavelength  components  of  actual  sea 
surface. 
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FIGURE  13.  Comparison  of  estimated  and  actual  time  histories  at  the  observation  locations. 

6  Conclusions  and  recommendations 

This  report  describes  a  variational  approach  for  assimilation  of  single-point  in-si tu  observation 
of  surface  elevation  and  velocity  into  the  extended  Boussinesq  model  of  Wei  et  al.  (1996).  The 
mathematical  framework  for  assimilation  is  presented,  along  with  applications  of  the  procedure  for 
two  cases:  monochromatic  waves  and  broad-spectrum  irregular  waves.  Overall  the  performance 
of  the  approach  for  monochromatic  waves  is  quite  good;  while  for  irregular  waves,  some  of  the 
short-wave  detail,  especially  for  waves  propagating  in  the  along-shore  direction,  is  missed.  As 
part  of  the  development  work  for  the  assimilation  procedure  a  parallel  Boussinesq  model  code 
was  improved  through  addition  of  run-up,  wave-breaking  and  and  bottom-friction  models,  and  a 
directional  wavemaker  was  implemented. 

The  assimilation  methodology  is  based  on  a  variational  approach  wherein  the  wavemaker  source 
function  for  the  Boussinesq  model  which  results  in  a  best  fit  between  the  calculated  wave  field 
and  the  observations  is  determined.  To  do  this,  a  cost  function  is  defined  which  is  a  positive- 
definite  measure  of  the  error  between  the  predictions  and  observations.  The  control  variable,  that 
which  is  adjusted  to  produce  the  desired  wave  field,  is  the  wavemaker  source  function  for  the 
Boussinesq  model.  To  find  the  minimum  in  the  cost  function,  subject  to  the  constraint  that  the 
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estimated  wave  field  satisfy  the  extended  Boussinesq  equations,  the  cost  function  is  augmented 
with  the  Boussinesq  equations  weighted  with  the  adjoint  variables.  Setting  the  first  variation  of  the 
augmented  cost  function  with  respect  to  the  solution  variables  and  the  adjoint  variables  results  in 
the  conditions  required  for  a  minimum;  it  also  produces  the  adjoint  equations.  The  solution  of  the 
adjoint  equations  is  used  to  calculate  the  gradient  of  the  cost  function  with  respect  to  the  control 
variable.  This  gradient  is  used  in  a  conjugate-gradient  iteration  scheme  to  find  the  wavemaker 
source  function  which  yields  the  best  fit  to  the  data  as  measured  by  the  cost  function.  Additional 
terms  related  to  the  wavemaker  source  function  are  included  in  the  cost  function  to  ensure  a  unique 
result. 

As  paxt  of  the  code  development  for  this  effort,  and  existing  parallel  Boussinesq  code  based 
on  both  the  formulation  and  numerical  approach  of  Wei  et  al.  (1996),  was  extended  to  include 
modeling  for  run-up,  wave-breaking  and  bottom  friction.  In  addition  both  monochromatic  and 
broad-spectrum  wavemakers  were  implemented.  This  results  in  a  Boussinesq  model  code  capable 
of  generating  realistic  wave  fields  either  as  a  product  of  the  assimilation  process,  or  as  a  direct 
forward  prediction  based  on  a  known  incident  wave  spectrum. 

The  assimilation  procedure  was  applied  to  two  cases,  observations  of  monochromatic  waves 
and  irregular  waves.  Application  of  the  procedure  for  monochromatic  waves  results  in  a  converged 
solution  after  about  five  conjugate  gradient  iteration  cycles.  Comparison  of  the  spatial  variations 
in  the  wave  field  at  an  instant  in  time  reveals  some  spatial  undulations  in  the  estimated  wave  field, 
but  otherwise  good  agreement.  This  indicates  that  some  spurious  higher- wave- number  components 
have  been  introduced  by  the  assimilation  procedure.  Comparison  of  the  estimated  time  series  at 
the  observation  locations  to  the  observations  shows  good  agreement  between  the  predictions  and 
observations. 

Use  of  the  assimilation  procedure  for  data  derived  from  broad-spectrum  irregular  waves  produces 
encouraging,  but  less-good  results.  The  assimilation  procedure  converges  in  about  ten  iterations, 
although  it  continues  to  drop,  albeit  slowly,  for  about  ten  more  iterations.  The  results  estimate  of 
the  wave  field  reproduces  the  longer  wavelength  components,  while  shorter  wavelength  components 
are  missed.  This  is  especially  true  for  the  along-shore  propagating  wave  components. 

The  approach  developed  and  presented  here  appears  promising,  allowing  the  entire  wave  field 
to  be  estimated  from  a  set  of  single-point  in-situ  observations.  The  resulting  estimate  is  consistent 
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with  the  physics  of  near-shore  waves  as  represented  in  extended  Boussinesq  models.  The  work 
under  this  effort  is,  however,  primarily  a  demonstration  that  such  an  approach  to  the  problem  is 
feasible.  Additional  necessary  work  includes: 

•  Exercising  the  assimilation  algorithm  for  a  broad  range  of  simulated  wave  fields  to  determine 
the  limitations  of  the  resulting  wave  estimates. 

•  Application  of  the  assimilation  algorithm  to  field  observation  and  comparison  to  other  ground- 
truth  data. 

•  Investigation  of  alternative  approaches  to  ensuring  uniqueness  of  the  resulting  wave  field.  The 
present  approach  penalizes  spatial  gradients  in  the  wavemaker  source  function  and  may  result 
in  suppression  of  along-shore  waves;  and  an  alternative  approach  which  does  not  have  this 
effect  would  be  desirable. 

Finally,  the  assimilation  algorithm  described  here  allows  one  to  relate  a  set  of  single-point  in- 
situ  observations  to  a  synoptic  picture  of  the  wave  field.  This  ability  can  most  likely  be  exploited 
to  develop  methods  to  determine  an  optimal  approach  for  observing  waves  in  a  give  region,  i.e. 
an  approach  which  extracts  the  maximum  amount  of  information  from  the  minimum  number  of 
observations. 
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